RELIABILITY ENHANCEMENT OF NAVIER-STOKES CODES 
THROUGH CONVERGENCE ENHANCEMENT 


K.-Y. Choi and G. S. Dulikravich 

Department of Aerospace Engineering, The Pennsylvania State University 
University Park, PA 16802 


SUMMARY 

Reduction of total computing time required by an iterative algorithm for solving Navier-Stokes equations is an 
important aspect of making the existing and future analysis codes more cost effective. Several attempts have been 

made to accelerate the convergence of an explicit Runge-Kutta time-stepping algorithm. These acceleration methods 
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are based on local time stepping , implicit residual smoothing , enthalpy damping and multigrid techniques 

Also, an extrapolation procedure based on the power method and the Minimal Residual Method (MRM) were 

applied to the Jameson’s multigrid algorithm. The MRM uses same values of optimal weights for the corrections to 

every equation in a system and has not been shown to accelerate the scheme without mul tigriding. Our Distributed 

Minimal Residual (DMR) method^ based on our General Nonlinear Minimal Residual (GNLMR) method^ allows 

each component of the solution vector in a system of equations to have its own convergence speed. The DMR 

method was found capable of reducing the computation time by 10-75% depending on the test case and grid used. 

Recently, we have developed and tested a new method termed Sensitivity Based DMR or SBMR method that is 

easier to implement in different codes and is even more robust and computationally efficient than our DMR method. 

TECHNICAL DISCUSSION 

This method predicts an optimum amount of correction to the solution vector by combining the information 
from several previous iteration levels. Each of the corrections obtained from the past iterations is multiplied by a 
different weighting factor and these weighting factdrs are determined so that they minimize the overall future 
residual. Although it is based on general Krylov subspace methods, the DMR method differs from them by the fact 
that weighting factors are different from one variable to another in the system. Recently, we divided the 
computational domain into several zones and the DMR method was applied separately in each of these zones. This 
approach did not show noticeable improvement over the original DMR method. Also it was difficult to decide how 
to divide the domain systematically into several zones. The DMR concept was also applied to minimize the future 
residual at each grid point rather than to minimize the residual integrated over the whole domain. This approach was 
not successful, because the weighting factors obtained for each grid point differ too much from one point to another 
thus making convergence history erratic and often diverging. 



SENSITIVITY BASED MINIMAL RESIDUAL (SBMR) METHOD 

The residual at a grid point depends on the solution vector Q at the neighboring points including the point 
itself. The sensitivity of residual R m (m=l,23: number of equations) with respect to Q k (k=typical neighboring 


points) is Si. Notice that the sensitivity can be determined from the finite difference equation used in the scheme 

aq 

For a two-dimensional incompressible flow solved using Chorin's artificial compressibility method the solution 
T 

vector is Q = [p u v] . Suppose we have calculated the solution vector Q at iteration levels up to t+n where n is 
the number of regular iteration steps between two iteration levels. Then the change in the solutions between the 
iteration levels can be written as follows. 

Ap k = (p k r-(pJ Au k = (u k r-(u J, Av k = (v k r-(vJ, (I) 

Using the first two terms in a Taylor series expansion in artificial time direction, each residual for a two-dimensional 
system after n iterations will be 
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Similarly, future residual at t = (t+n)+l can be approximated by 
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Here, a's are the factors that multiply A's to estimate the future solution vector so that it satisfy a desired objective. 
For now, each a is assumed to have the same value over the whole domain. The future solutions are estimated as 

(p/'^IpJf+aAp, (4) 

with similar expressions for u and v. The a’s are determined such that the L-2 norm of the overall future residual will 
be minimized, that is. 
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where the subscript q stands for each flow variable p, u and v. For simplicity, let us denote the bracketed terms in 
equation (3) as a , a and a so that 
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with similar expressions for a and a . Therefore, 
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Substituting (7) into (5) gives the following three equations for optimal global a’s. 
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with similar expressions involving and U mv . In equation (8), R's and a’s are known from the past iteration 

results. Since each a is assumed to have the same value over the computational domain, equation (8) gives a system 

of simultaneous equations for a , a and a . 
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with similar expressions arising from u and v components. As the grids are clustered (higher cell aspect ratio), local 
time steps become smaller in those regions in order to meet the stability criterion. Numerical results (Figures 1-6) 
obtained at low and high Reynolds numbers on non-clustered and moderately clustered grids for straight and U- 
shaped two-dimensional channel flow demonstrate that SBMR method can: a) be used in conjunction with any basic 
iterative algorithm, b) be used with only minor modifications in the existing codes, c) significantly accelerates 
iteration procedure, d) perform more effectively at lower Reynolds numbers, e) perform well on moderately clustered 
grids. 

The SBMR method applied so far calculates the same a f s for the entire computational domain, which cannot 
represent optimum a's for both coarse grid regions and fine grid regions. Therefore, a natural conclusion is to allow 
a's to have different values in the clustered regions. The future work on SBMR will concentrate on applying this 
method by lines such that each grid line normal to the clustering direction has its own a's. We expect that with this 
approach, SBMR concept will calculate the optimal local corrections to solution vectors on highly clustered grids. 
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Figure 1. Convergence histories for a cascade flow: Re=500, non- 
clustered 40x30 grid, CFL=2.8, p=3.0. No artificial dissipation was 
used Half of NACA 0012 airfoil on top and bottom walls. 



Figure 2. Convergence history for a straight channel flow: Re=1600, non- 
clustered 60x60 grid, x-length=5, y-length=l, p=5, von=0.4, CFL=2.5, 
omega=0.001. Initial guess: u=1.0E-5, v=0, p=0. Inlet velocity profile: 
u= parabolic, v=0. 
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Figure 3. Convergence history for a straight channel flow: Re=1.6 million, non- 
clustered 60x60 grid, x-length=5, y-Iength=l, p=5, CFL=2.8, von=0.4“^, 
omega=0.005. Initial guess: u=1.0E-5, v=0, p=0. Inlet velocity profile: 
u=parabolic, v=0. 



Figure 4. Straight channel flow: Re=1600, max AR=100 on 60x60 grid 
x-length=5, y-length=l, 0=5, CFL=2.5, von=0.4, omega=0.0. Initial guess: 
u=1.0E-5, v=0, p=0 everywhere. Inlet velocity profile: u=parabolic, v=0. 
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Figure 5. Straight channel flow: Re=1.6 million, max AR=100 on 60x60 grid, 

x-length=5, y-length=l, p=5 t CFL=2.5, von=0.4xl0^, omega=0.005. Initial 
guess: u=1.0E-5, v=0, p=0. Inlet velocity profile: u=parabolic, v=0. 



Figure 6. Convergence history for a U-shaped channel flow (Re=100, 
129x30 grid cells); f£=5, CFL=2.8, von=0.4, omega=0.0. 


I 


186 










